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INTRODUCTION 


Our  originsl  proposal  was  in  two  parts.  Tho  first  involved  improved  optimization 
methods  for  clusters  and  weakly  interacting  systems.  Traditional  optimization  methods 
do  not  take  into  account  the  special  character  of  these  systems,  such  as  the  large 
disparity  between  the  forces  within  and  between  molecules,  and  are  thus  inefficient. 
Our  aim  here  was  to  develop  an  efficient  optimization  procedure  for  weakly  interacting 
molecular  clusters,  as  well  as  for  other  non-standard  optimization  scenarios,  such  as 
molecular  adsorption  on  a  model  surface. 

The  second  project  dealt  with  improved  electronic  structure  methods  for  large  molecular 
systems.  The  main  thrust  here  was  the  development  of  efficient  algorithms  for  the 
computation  of  both  canonical  and  local  MP2  energies  and  gradients.  Current  density 
functionals  are  not  able  to  adequately  model  dispersion  forces  and  MP2  should  be  a 
good  option  here.  These  latter  calculations  can  eventually  supply  molecular  energies 
and  gradients  for  the  optimization  project.  We  have  also  developed  a  preliminary 
version  of  an  atom-centered  fast  multipole  program  for  large-molecule  DFT 
calculations. 

This  final  report  summarizes  the  progress  made  in  these  areas  during  the  period  of  this 
grant  award.  No  great  detail  is  given,  as  this  has  to  a  large  extent  already  been 
provided  in  our  previous  yearly  reports. 


1 .  Optimization  of  Molecular  Geometries 
a.  Cluster  and  Rigid-Body  Optimization 

We  have  successfully  developed  and  fully  implemented  an  algorithm  for  the  efficient 
optimization  of  molecular  clusters.  This  uses  a  special  set  of  what  we  have  called 
cluster  coordinates,  which  comprise  the  usual  stretches,  bends  and  torsions  to  describe 
the  interactions  within  each  molecule  (the  intramolecular  geometry)  and  inverse- 
distance  coordinates  to  describe  the  interactions  between  individual  molecules  in  the 
cluster  (the  intermolecular  geometry).  Our  new  cluster  coordinates  reduce  the  number 
of  geometry  optimization  cycles  required  for  convergence  by  up  to  an  order  of 
magnitude  relative  to  Cartesian  coordinates,  and  also  significantly  relative  to  other 
coordinates,  e.g.,  distance  coordinates.  Molecular  clusters  usually  have  a  number  of 
minima,  and  cluster  coordinates  are  often  more  successful  in  reaching  low-lying  minima 
than  are  Cartesians. 

In  addition,  by  using  a  Schmidt-Orthogonalization  scheme  for  imposing  constraints 
introduced  earlier  (J.Baker,  A.Kessi  and  B.Delley,  J.Chem.Phys.  105  (1996)  192),  it  is 
possible  to  effectively  constrain  all  the  intramolecular  degrees  of  freedom  in  a  molecular 
cluster  and  in  this  way  carry  out  complete  rigid-body  optimizations. 

We  illustrate  the  efficacy  of  our  new  cluster  coordinates  by  optimizing  the  geometries  of 
twenty  randomly-generated  clusters,  each  containing  ten  hydrogen  molecules.  Table  1 
shows  the  starting  energy,  the  number  of  optimization  cycles  needed  to  converge  and 
the  final  energy  for  optimizations  using  both  Cartesian  and  cluster  coordinates.  As  can 
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be  seen,  the  average  rate  of  convergence  is  over  ten  times  faster  with  cluster 
coordinates  than  with  Cartesians,  and  in  each  case  the  optimization  converged  to  a 
lower  final  energy. 


Starting  energy,  number  of  cycles  to  converge  and  final  energy  for  the  RHF/3-21G 
optimization  of  20  randomly  generated  clusters  of  10  H2  molecules  (initial  bond  length 
0.72  A) 


cluster 


1 

starting  energy 

-11.166613 

2 

-11.151158 

3 

-11.159068 

4 

-11.142043 

5 

-11.177657 

6 

-11.178220 

7 

-11.150748 

8 

-11.152618 

9 

-11.170617 

10 

-11.148964 

11 

-11.157261 

12 

-11.171696 

13 

-11.161546 

14 

-11.164816 

15 

-11.143458 

16 

-11.151136 

17 

-11.157329 

18 

-11.171771 

19 

-11.141685 

20 

-11.150264 

average 


Cartesian 


cycies 

energy^* 

613 

-11.230147 

506 

-11.230130 

282 

-11.230023 

100 

-11.229742 

776 

-11.230133 

305 

-11.229987 

338 

-11.230103 

410 

-11.230122 

454 

-11.230144 

425 

-11.230093 

431 

-1 1 .230073 

501 

-11.230054 

500 

-11.230130 

435 

-11.230041 

478 

-11.230110 

965 

-11.230106 

555 

-11.230116 

357 

-11.230074 

172 

-11.229711 

307 

-11.230016 

446 

-11.230053 

cluster  coordinates^ 


cycles 

energy** 

34 

-11.230188 

33 

-11.230208 

30 

-11.230193 

55 

-11.230162 

24 

-11.230170 

58 

-11.230176 

31 

-11.230226 

47 

-11.230212 

17 

-11.230226 

42 

-11.230196 

70 

-11.230190 

39 

-11.230223 

66 

-11.230223 

33 

-11.230176 

32 

-11.230226 

54 

-11.230208 

33 

-11.230206 

41 

-11.230190 

39 

-11.230206 

35 

-11.230202 

41 

-11.230200 

®  intermolecular  cluster  coordinates  were  based  on  simple  inverse  distances  (1/R) 

^  convergence  criteria  were  a  maximum  gradient  component  <  0.00005  au  and  an  energy 
change  of  <  10'^  hartree 
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b.  Geometry  Optimization  for  Reactions  on  Modei  Surfaces 

We  have  developed  an  efficient  and  automatic  procedure  for  optimizing  the  geometries 
of  surface-adsorbate  systems.  Internal  coordinates  (stretches,  bends  and  torsions)  are 
generated  for  the  surface  and  the  molecule  being  adsorbed  separately,  and  additional 
linking  primitives  are  generated  to  connect  the  two  parts.  The  often  high  atomic 
connectivity  in  the  model  surface  makes  it  possible  to  ignore  torsions,  spanning  all 
surface  degrees  of  freedom  using  stretches  and  bends  only,  greatly  reducing  the 
dimension  of  the  primitive  space.  Cartesian  coordinates  (i.e.,  fixed  atoms)  can  be 
constrained  within  a  delocalized  internal  coordinate  optimization  (by  constraining  all 
primitives  that  involve  solely  the  fixed  atoms  in  their  definition),  and  the  effects  of 
freezing  various  surface  layers  on  reaction  energetics  and  geometries  can  be 
investigated. 

We  have  used  this  procedure  to  investigate  the  dissociative  adsorption  of  water  on  the 
reactive  site  of  the  silicon  (100)  surface,  modelled  using  Si9Hi2  and  SiirHao  clusters 
(see  Fig.  1).  The  Si  atoms  labelled  1  and  2  in  both  cases  constitute  the  "active  site"  and 
these  two  atoms  are  formally  unsaturated.  They  are  considered  to  be  in  the  first  (or 
surface)  layer  in  our  model.  Si  atoms  3,  4,  5  and  6  comprise  the  second  layer.  All  other 
Si  atoms  are  taken  as  the  third  layer.  Apart  from  the  two  surface  atoms,  the  remaining 
Si  atoms  are  saturated  by  adding  hydrogens. 

The  reaction  is  considered  to  take  place  by  first  forming  a  complex  between  water  and 
the  Si  surface,  with  a  bond  forming  between  the  water  oxygen  and,  say,  Si  atom  1 ,  with 
Si  atoms  1  and  2  and  the  water  molecule  all  in  the  same  plane;  this  is  followed  by  a 
lengthening  of  one  of  the  water  0-H  bonds  and  the  partial  formation  of  a  Si-H  bond  with 
Si  atom  2  (in  the  transition  state),  leading  to  the  final  product  with  the  water  molecule 
completely  dissociated  and  full  Si-OH  and  Si-H  bonds  at  the  previously  unsaturated  Si 
centers.This  process  is  shown  schematically  in  Fig.  2.  By  systematically  fixing  Si  atoms 
(and  their  attached  hydrogens)  in  the  three  surface  "layers"  defined  above,  we  can 
investigate  the  effects  on  the  reaction  energetics  of  freezing  various  surface  layers. 

Reaction  energetics  are  shown  in  Table  2.  We  used  density  functional  theory 
(specifically  the  hybrid  B3LYP  functional  with  the  standard  6-31 G*  basis  set).  All 
structures  were  taken  to  be  closed  shell  singlets.  The  initial  geometries  of  both  Si9Hi2 
and  Sii7H2o  were  fully  optimized  at  B3LYP/6-31G*. 
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Table  2  Energetics  of  the  dissociative  adsorption  reaction  (H2O  — >  H  +  OH)  on  the 
model  Si9Hi2  and  Sii7H2o  surfaces  (in  kcal/mol  relative  to  the  initial  adsorbed  complex 
as  the  energy  zero) 


surface  type 
fixed  surface 
layer  1  free 
layers  1  &  2  free 
unconstrained 

surface  type 

fixed  surface 
layer  1  free 
layers  1  &  2  free 
unconstrained 


Si9Hi2  +  H2O 
+5.3 
+11.7 
+14.4 
+14.9 

H2O-SI9H12 

0.0 

0.0 

0.0 

0.0 

transition  state 
+3.1 
+4.0 
+5.3 
+5.5 

SigHisOH 

-53.0 

-53.8 

-52.3 

-53.0 

SI17H20  ^  H2O 
+5.8 
+12.9 
+15.1 
+15.8 

H2O-SI17H20 

0.0 

0.0 

0.0 

0.0 

transition  state 
+2.8 
+4.0 
+5.1 
+5.1 

Sii7H2iOH 

-53.1 

-53.8 

-52.4 

-52.7 

As  can  be  seen  from  Table  2,  the  effect  on  the  overall  relative  energetics  of  freezing 
various  surface  layers  is  minor.  The  initial  binding  between  water  and  the  model  surface 
is  fairly  strong  systematically  increases  with  increasing  surface  relaxation,  this  must  be 
the  case  as  we  are  comparing  fully  optimized  reactants  with  only  a  partially  optimized 
complex  -  relaxation  of  each  surface  layer  can  only  lower  the  energy  of  the  complex, 
thus  increasing  the  binding.  However,  the  relative  energy  difference  between  the 
adsorbed  complex,  the  transition  state,  and  the  dissociatively  adsorbed  product,  hardly 
differs  as  we  systematically  relax  the  surface.  The  only  real  trend  is  a  small  but  steady 
increase  in  the  transition  state  barrier  with  increasing  surface  relaxation.  The  energetics 
on  the  totally  frozen  compared  to  the  totally  relaxed  surface  are  very  similar,  as  are  the 
results  for  the  Sii7H2o  cluster  compared  to  Si9Hi2. 

Given  the  quality  of  the  calculations  (a  reasonable  but  not  especially  high  level  of 
theory),  our  conclusions  are:  (1)  that  the  frozen  surface  model  provides  a  pretty  good 
picture  of  the  overall  reaction  energetics;  and  (2)  that  the  Si9Hi2  model  surface  seems 
to  be  perfectly  adequate  to  describe  the  dissociative  adsorption  of  a  water  molecule  on 
a  silicon  cluster. 
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Figure  2:  Schematic  of  geometrical  arrangement  around  the  active  site  for  the 
dissociative  adsorption  of  water  on  a  model  Si  (100)  surface  for:  (a)  the  initial  adsorbed 
complex:  (b)  the  transition  state;  (c)  the  dissociatively  adsorbed  product  (see  Table  1) 
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c.  Large-Molecule  Optimization 

The  bottleneck  preventing  the  efficient  use  of  internal  coordinates  in  large  systems  is 
the  time  taken  to  transform  the  gradient  (and  possibly  the  Hessian)  from  Cartesian 
coordinates  to  internal  coordinates,  and  to  transform  the  new  geometry  from  internal 
coordinates  back  to  Cartesians.  The  latter  step  is  needed  to  calculate  the  energy  and 
gradient  for  the  next  optimization  step.  Both  these  transformation  steps  formally  scale 
as  the  cube  of  the  system  size,  O(N^).  We  have  developed  several  novel  algorithms  to 
reduce  this  scaling  to  near-linear. 

The  gradient  transformation  can  formally  be  written  as 

gCart  _  gtgint  (1) 

where  B  is  the  (linearized)  transformation  matrix  from  Cartesian  to  internal  coordinates 
(the  so-called  Wilson  B-matrix).  As  it  is  not  a  square  matrix,  does  not  have  a  regular 
inverse,  and  the  usual  trick  is  to  define  (B^)  ^  =  (BB^)  B,  which  acts  as  a  left  inverse  in 
the  sense  that  (bVb^=  1.  The  cubic  computational  dependence  of  traditional  internal 
coordinate  transformations  is  caused  by  the  presence  of  this  formal  matrix  inverse. 


Using  this  left  inverse  we  can  write 

gint  _  =  (BBVBq^'^  (2) 

which  can  be  rearranged  to 

(BB^)g‘"*  =  Bg“^  (3) 

Splitting  off  the  diagonal  D  of  BB^,  we  can  rewrite  eq.  3  as 

D  g'"‘  =  B  g*^"*  +  (D  -  BB^)  g‘"*  (4) 

and  convert  it  into  an  iterative  scheme 

g"’t(k+l )  =  D-'  [  B  g“^  +  (D  -  BB^)  g‘"‘(k)  ]  (5) 


which  can  be  solved  by  the  preconditioned  conjugate  gradient  method.  Solution  of  Eq. 
5  is  formally  O(N^)  as  it  no  longer  contains  a  matrix  inverse.  In  practice,  however,  the 
actual  order  depends  on  how  rapidly  the  procedure  converges:  if  it  takes  0(N)  cycles  to 
converge,  then  nothing  is  gained. 

The  backtransformation  of  the  new  geometry  from  internal  to  Cartesian  space  is  also 
accomplished  iteratively  via 

x(k+1)  =  x(k)  +  (BV(k)[q  ■  q(k)]  (6) 

Here  q  is  the  (known)  new  geometry  in  internal  coordinates  and  the  iterative  procedure 
is  typically  started  (k=1)  using  the  old  Cartesian  geometry  (and  the  old  B  matrix). 
Although  this  iterative  back  transformation  normally  converges  extremely  rapidly 
(differences  between  q  and  q(k)  of  <  10'^°  in  3-5  cycles),  it  formally  requires  the 
construction  and  inversion  of  a  new  B  matrix  at  each  iterative  cycle.  A  better  scheme, 
which  works  unless  the  displacements  are  very  large,  is  to  keep  the  old  (B^)'  matrix 
throughout;  this  increases  the  number  of  cycles  but  each  cycle  is  now  much  less 
expensive  as  no  matrix  inverse  is  needed. 
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Our  first  approach  involved  direct  solution  of  Eq.  2  using  natural  internal  coordinates 
and  a  scheme  in  which  the  formal  inverse  of  BB^  is  replaced  by  the  inverse  of  its 
approximate  symmetric  Cholesky  decomposition  LdL^,  where  L  is  lower  triangle  and  d 
is  diagonal.  In  the  decomposition,  we  keep  only  elements  of  L  that  exceed  a  small 
threshold;  this  number  is  expected  to  grow  approximately  linearly  with  system  size.  This 
reduces  the  gradient  transformation  to  roughly  O(N^);  it  also  has  the  advantage  that  the 
same  decomposition  can  be  used  for  the  backtransformation,  Eq.  6  (assuming  that  the 
geometry  change  is  not  too  large). 

We  subsequently  developed  an  alternative  approach  using  delocalipd  internal 
coordinates  and  the  iterative  gradient  transformation  of  Eq.  5.  Delocalized  internals  are 
formed  by  diagonalization  (or  its  equivalent)  of  BB^  on  the  first  optimization  cycle  and, 
aftGr  transformation  to  tho  new  coordinate  space,  the  BB^  matrix  is  usually  diagonally 
dominant  throughout  the  optimization.  This  ensures  that  Eq.  5  converges  rapidly. 
However,  the  main  advantage  of  delocalized  internals  is  that  the  formal  0(N  )  traditional 
backtransformation  can  be  replaced  by  a  Z-matrix  approach  which  is  0(N).  (This  Z- 
matrix  backtransformation  cannot  readily  be  used  with  natural  internals  because  they 
mix  in  too  much  of  the  formally  redundant  part  of  the  full  space  of  primitive  internals.) 
Unfortunately,  delocalized  internals  -  being  linear  combinations  of  all  possible  stretches, 
bends  and  torsions  in  the  primitive  space  (unlike  natural  internals,  which  are  linear 
combinations  of  just  a  few  primitives  and  are  highly  localized)  -  are  not  efficient  beyond 
a  certain  molecular  size  (around  500  or  so  atoms)  due  to  their  lack  of  sparsity.  Natural 
internals  are  sparse,  but  are  increasingly  difficult  to  generate  as  molecular  size  and 
topological  complexity  increases. 

Our  most  recent  approach  is  capable  of  carrying  out  large-molecule  geometry 
optimizations  directly  in  the  full  space  of  primitive  internal  coordinates.  Normally  far 
more  primitive  stretches,  bends  and  torsions  can  be  generated  for  a  large  molecule 
containing  N  atoms  than  the  3N-6  that  are  necessary  to  span  all  the  internal  degrees  of 
freedom.  Although  the  B  matrix  in  primitive  internals  can  have  a  much  greater 
dimension  than,  say,  its  equivalent  in  delocalized  internals,  it  is  highly  sparse,  having  at 
most  12  non-zero  entries  per  column  regardless  of  the  number  of  primitives,  which 
could  easily  be  in  the  tens  of  thousands  for  a  system  with  a  few  thousand  atoms.  If  we 
let  g'"*  =  Bx,  then  substituting  into  Eq.  1  gives 

gCart  _  gtgx  (7) 

We  can  solve  this  set  of  linear  equations  to  get  x,  and  then  back-substitute  to  get  g'"*. 
This  involves  the  formal  inverse  of  B^B  (instead  of  BB^),  and  not  only  does  this  matrix 
have  smaller  dimension  than  BB^  (3N-6  instead  of  the  size  of  the  primitive  space),  but  it 
turns  out  that  its  approximate  Choleskly  decomposition  has  a  sparser  L  matrix  than  for 
the  corresponding  BB'*'  matrix,  leading  to  a  further  reduction  in  CPU  time  for  the 
transformation. 

Work  in  this  area  is  ongoing,  generously  supported  by  a  new  Air  Force  grant  (Award 
number  F49620-00- 1-0281),  and  we  hope  to  further  improve  the  performance  of  our 
large  molecule  optimization  algorithm. 

To  give  some  idea  of  the  efficiency  of  our  coordinate  transformations,  and  the 
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performance  of  the  optimization  overall,  we  present  the  following  three  Tables.  Table  3 
shows  timings  for  various  steps  in  the  large-molecule,  delocalized  internal  coordinate 
algorithm,  together  with  the  timings  (old)  if  all  transformations  were  done  using  full 
matrix  inversions. 


Timings  (CPU;  seconds)  on  an  IBM  RS6000/390  workstation  for  various  steps  in  the 
internal-coordinate  large-molecule  optimization  algorithm. 


cycle®  B  Matrix 

construction 


gradient^ 

transformation 


back® 

transformation 


old 

new 

old 

new 

old 

new 

Jawsamycin  (C32H43N3O6,  84 

atoms) 

16.54 

0.12 

1  34.01 

12.61 

2.01 

0.03 

2  2.09 

0.27 

2.03 

0.10 

16.54 

0.14 

3  2.12 

0.27 

2.02 

0.12 

16.54 

0.17 

30  2.12 

0.27 

2.01 

0.55 

12.37 

0.23 

Taxol  (C47H51NO14 

,113  atoms) 

36.19 

0.23 

1  110.06 

23.09 

4.46 

0.05 

2  4.57 

1.45 

4.41 

0.22 

36.15 

0.25 

3  4.58 

1.45 

4.41 

0.22 

36.15 

0.25 

68  4.59 

1.44 

4.46 

1.36 

27.10 

0.98 

(Alanine)2o  (C60H102N20O21,  203  atoms) 

492.97 

0.42 

1  826.07 

97.25 

99.91 

0.29 

2  23.71 

9.30 

99.96 

1.07 

626.28 

0.51 

3  23.80 

9.27 

100.01 

1.42 

620.58 

0.51 

63  23.80 

9.29 

100.37 

2.50 

372.89 

0.58 

total 

time 

old  new 


53.41  13.61 

20.68  0.53 

20.70  0.58 

16.52  1.07 


152.24  24.90 

45.17  1.96 

45.18  1.96 

36.19  3.82 


1439.47  118.48 
750.13  11.06 

744.57  11.38 

497.24  12.55 


®  Shown  are  timings  for  the  first  three  optimization  cycles  and  the  final  cycle.  Convergence  criteria 

are;  maximum  gradient  component  <  0.00005  Ef/aoi  energy  change  <  10  Eh 

®  Timings  for  B  matrix  construction  on  cycle  1  include  once  only  construction  of  delocalized  internals 
The  fast  gradient  and  back  transformation  were  converged  to  an  accuracy  of  better  than  5x10 
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As  can  be  seen,  the  savings  are  enormous  even  for  these  relatively  small  systems.  Our 
most  recent  B^B  algorithm,  using  either  natural  internal  coordinates  or  a  full  set  of 
primitive  internals,  is  even  faster  for  all  steps  except  possibly  the  backtransformation. 
Table  4  shows  timings  for  the  gradient  transformation  and  the  total  (average)  time  per 
optimization  cycle  for  a  series  of  alanine  polypeptides;  also  shown  are  the  density  (the 
percentage  of  non-zero  elements)  in  the  Cholesky  decomposition  for  both  the  B  B  and 
the  original  BB^  algorithm. 


Table  4.  Comparison  of  the  performance  of  the  B^B  and  BB^  formalisms  for  alpha- 
helical  alanine  polypeptides.  The  following  data  are  shown:  the  density  (the  fraction  of 
non-zero  elements)  in  the  Cholesky  factors,  and  the  average  time  (Tcydei  sec)  per 
optimization  cycle  (a  few  decompositions,  one  force  transformation,  and  a  line  search 
including  4-5  geometry  back-transformations)  on  a  300  MHz  Pentium  II  PC  for  both  the 
B'^B  and  BB''’  formalisms.  For  the  B^B  formalism  the  average  time  (Tporce,  sec)  of  a 
single  force  transformation  step  is  also  shown. 


Alanine  polypeptides 

#  of  alanine 

#  of  atoms 

Density  % 

Tporce 

Tcycle 

Density  % 

Tcycle 

units 

10 

109 

17.95 

0.05 

0.77 

26.37 

1.09 

20 

209 

9.63 

0.11 

1.61 

18.86 

2.70 

30 

309 

6.58 

0.21 

2.62 

13.20 

5.21 

40 

409 

4.99 

0.24 

3.88 

10.77 

8.24 

50 

509 

4.02 

0.37 

5.23 

9.32 

12.29 

100 

1009 

2.04 

0.70 

11.17 

5.94 

40.04 

200 

2009 

1.03 

1.33 

22.74 

3.88 

137.97 

300 

3009 

0.69 

2.01 

.  32.28 

- 

500 

5009 

0.41 

2.49 

49.55 

- 

- 

700 

7009 

0.30 

3.40 

71.45 

- 

- 

999 

9999 

0.21 

4.57 

91.62 

- 

- 

Finally  Table  5  directly  compares  the  performance  of  our  large-rnolecule  internal 
coordinate  optimization  algorithm  with  the  standard  conjugate  gradient  optimizer  in 
Cartesian  coordinates  in  the  TINKER  molecular  mechanics  package  (J.W. Ponder, 
Software  Tools  for  Molecular  Design,  version  3.7,  June  1999)  for  the  same  alanine 
polypeptides  shown  in  Table  4.  As  can  be  seen,  our  algorithm  converges  in  dramatically 
fewer  energy  and  gradient  evaluations  than  the  more  traditional  Cartesian  optimizer  and 
typically  converges  to  a  lower  final  energy.  This  is  particularly  apparent  for  the  larger 
systems. 
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Table  5  Comparison  of  optimizations  of  alpha-helical  alanine  polypeptides  using 
natural  internal  and  Cartesian  coordinates.  The  initial  and  final  energies  (kcal/rnol)  and 
the  number  of  gradient  and  energy  evaluations  are  shown  for  Cartesian  optimizations 
as  implemented  in  the  TINKER  program  {minimize.x)  and  natural  internal  coordinate 
optimizations  in  our  program 


Natural  Internal 


#of 

initial 

#of 

alanines 

energy 

E 

G 

10 

11.4867 

31 

138 

20 

6.2427 

38 

168 

30 

0.7077 

50 

216 

40 

-4.8271 

53 

237 

50 

-10.3619 

64 

314 

100 

-38.0366 

71 

340 

200 

-427.7795 

69 

336 

300 

-650.1882 

76 

368 

500 

-2630.4851 

48 

234 

700 

-3688.1297 

50 

244 

999 

-5269.9224 

55 

267 

Coordinates 

TINKER  (0.01  kcal/mol  A) 

final 

#of 

final 

energy 

E 

G 

energy 

-43.3114 

313 

578 

-43.3124 

-97.4764 

412 

770 

-97.4768 

-152.1734 

354 

675 

-152.1717 

-206.8709 

335 

636 

-206.8696 

-261.5697 

436 

797 

-261.5688 

-535.0629 

664 

1246 

-535.0583 

-1082.0471 

1589 

3020 

-1081.9656 

-1629.0323 

2191 

4051 

-1628.7909 

-2723.0002 

6369 

12228 

-2719.7795 

-3816.9708 

6536 

12614 

-3799.8068 

-5452.4565 

7851 

15165 

-5414.0775 
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II.  Electronic  Structure  Theory 
a.  Local  Electron  Correlation 

We  have  continued  developing  our  local  second-order  Moller-Plesset  perturbation 
theory  code. 

<This  part  to  be  completed  by  PP> 


b.  Atom-centered  fast  multipoles 

The  aim  of  this  research  project  is  to  develop  a  fast  and  accurate  DPT  program  for 
large  molecules.  Instead  of  calculating  all  the  ^o-electron  integrals  as  in  the  traditional 
quantum  chemical  approach,  we  expand  the  Coulomb  potential  in  terms  of  spherical 
harmonics  over  atom-centered  grids.  Our  approach  is  similar  to  that  of  Delley  in  the 
DMol  program  (B.Delley,  J.Chem.Phys.  92  (1990)  508),  except  that  we  subsequently  fit 
the  potential  rather  than  do  a  global  fit  to  the  density  as  is  the  case  in  DMol.  In  large 
molecules,  the  Coulomb  potential  between  atoms  that  are  far  apart  can  be  described  by 
a  simple  multipole  expansion,  and  this  is  where  we  hope  to  gain  efficiency  over  the 
more  traditional  methods.  Our  approach  is  conceptually  similar  to  the  Continuous  Fast 
Multipole  Expansion  (C.A.White,  B.G. Johnson,  P.M.W.Gill  and  M.Head-Gordon, 
Chem.Phys.Lett.  230  (1994)  8).  However,  in  CFMM  the  molecule  is  divided  into  artificial 
rectangular  boxes,  with  multipole  expansions  being  used  to  describe  the  Coulomb 
interactions  between  boxes  that  are  far  enough  apart.  These  rectangular  boxes  require 
high  order  multipoles  and  are  responsible  for  the  fact  that  CFMM  becomes  really 
efficient  only  for  huge  systems.  Our  approach,  with  spherical  atom-centered 
expansions,  is  more  physical  and  expected  to  give  a  crossover  with  traditional  methods 
at  much  smaller  system  size. 

Initially  we  have  concentrated  on  producing  accurate  and  reliable  DFT  energies  for 
small  and  mediurn-sized  molecules  and  have,  we  believe,  obtained  some  of  the  most 
accurate  energies  using  this  approach  that  have  been  reported  to  date.  We  have  been 
very  careful  with  our  numerical  integration  techniques,  and  can  take  the  various 
expansions  involved  to  high  orders  in  spherical  harmonics  to  reach  limiting  values. 
Using  identical  Gaussian  basis  sets,  we  can  directly  compare  energies  from  the  fully 
numerical  spherical  harmonics  code  with  energies  from  traditional  quantum  chemistry 
DFT  codes  to  determine  how  high  to  take  L  in  the  spherical  harmonics  expansion  to 
accurately  reproduce  traditional  DFT  energies.  Our  findings  suggest  that  high  values  of 
L  are  needed  even  to  reproduce  relative  energies  accurately,  and  that  the  default 
expansions  in  many  commonly  used  codes  are  probably  too  low  to  give  reliable 
energetics. 

Table  6  shows  energies  for  a  number  of  di-  and  tri-atomics  for  L  values  in  the  spherical 
harmonics  expansion  from  0  to  19.  For  most  of  the  molecules  investigated,  absolute 
energies  do  not  stabilize  until  at  least  L=11  and  for  some,  e.g.  F2,  not  even  then. 
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Table  6  Analytical  and  numerical  BLYP/6-31G**  SCF  energies  (hartree)  for  a  number 
of  small  molecules  for  L  values  in  the  spherical  harmonic  expansions  from  0  to  19 


Lmax 

H2 

LiH 

0 

-1.159720 

-8.099536 

1 

-1.156331 

-8.070451 

2 

-1.166369 

-8.068231 

3 

-1.168114 

-8.064447 

4 

-1.168037 

-8.064419 

5 

-1.167918 

-8.066075 

6 

-1.167889 

-8.066663 

7 

-1.167888 

-8.066573 

8 

-1.167890 

-8.066523 

9 

-1.167890 

-8.066545 

10 

-1.167890 

-8.066544 

11 

-1.167890 

-8.066533 

12 

-1.167890 

-8.066531 

13 

-1.167890 

-8.066531 

14 

-1.167890 

-8.066530 

15 

-1.167890 

-8.066530 

16 

-1.167890 

-8.066530 

17 

-1.167890 

-8.066530 

18 

-1.167890 

-8.066530 

19 

-1.167890 

-8.066530 

no  pruning* 

-1.167891 

-8.066535 

exact 

-1.167896 

-8.066464 

HF 

Liz 

LiF 

-100.552673 

-15.158289 

-108.390760 

-100.321579 

-14.965681 

-107.743656 

-100.376894 

-14.979977 

-107.459285 

-100.409630 

-14.991192 

-107.380373 

-100.414410 

-14.994298 

-107.371944 

-100.412624 

-14.993896 

-107.382928 

-100.411409 

-14.993019 

-107.393007 

-100.411075 

-14.992559 

-107.398321 

-100.411025 

-14.992445 

-107.400443 

-100.411029 

-14.992465 

-107.401152 

-100.411050 

-14.992502 

-107.401342 

-100.411075 

-14.992524 

-107.401334 

-100.411078 

-14.992527 

-107.401307 

-100.411080 

-14.992527 

-107.401258 

-100.411079 

-14.992527 

-107.401204 

-100.411079 

-14.992527 

-107.401158 

-100.411079 

-14.992527 

-107.401126 

-100.411079 

-14.992527 

-107.401107 

-100.411079 

-14.992527 

-107.401100 

-100.411079 

-14.992527 

-107.401098 

-100.411094 

-14.992534 

-107.400661 

-100.411080 

-14.992538 

-107.400627 
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Table  6  (continued) 


Lmax 

CO 

N2 

Fz 

HzO 

CO2 

0 

-112.815305 

-108.970494 

-199.188346 

-76.444437 

-187.634789 

1 

-113.259613 

-109.387299 

-199.355946 

-76.249007 

-188.254279 

2 

-113.207833 

-109.412242 

-199.500291 

-76.365904 

-188.444913 

3 

-113.269420 

-109.490394 

-199.476625 

-76.404018 

-188,490216 

4 

-113.292217 

-109.512858 

-199.480261 

-76.403080 

-188.569540 

6 

-113.296366 

-109.515300 

-199.488192 

-76.399434 

-188.566628 

6 

-113.295616 

-109.513217 

-199.492640 

-76.398689 

-188.566551 

7 

-113.294325 

-109.511063 

-199.493785 

-76.398819 

-188.565544 

8 

-113.293430 

-109.509823 

-199.493646 

-76.398901 

-188.563028 

9 

-113.292972 

-109.509398 

-199.493387 

-76.398906 

-188.563087 

10 

-113.292799 

-109.509394 

-199.493261 

-76.398917 

-188.562201 

11 

-113.292770 

-109.509498 

-199.493201 

-76.398933 

-188.562271 

12 

-113.292775 

-109.509517 

-199.493181 

-76.398935 

-188.562270 

13 

-113.292783 

-109.509524 

-199.493160 

-76.398934 

-188.562268 

14 

-113.292787 

-109.509525 

-199.493146 

-76.398934 

-188.562287 

15 

-113.292789 

-109.509525 

-199.493139 

-76.398933 

-188.562287 

16 

-113.292790 

-109.509524 

-199.493138 

-76.398933 

-188.562290 

17 

-113.292790 

-109.509524 

-199.493139 

-76.398932 

-188.562291 

18 

-113.292789 

-109.509524 

-199.493142 

-76.398932 

-188.562290 

19 

-113.292789 

-109.509524 

-199.493143 

-76.398932 

-188.562290 

no  pruning* 

-113.292929 

-109.509526 

-199.492959 

-76.398934 

-188.562392 

exact 

•113.292888 

-109.509547 

-199.492937 

-76.398886 

-188.562330 

*  the  various  atom-centered  grids  are  normally  “pruned”  to  reduce  the  number  of  angular  grid  points  in 
appropriate  regions  of  space,  e.g.,  close  to  the  nucleus;  “no  pruning”  energies  are  with  the  full  grid 
without  any  reduction  in  the  number  of  angular  grid  points. 
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To  give  some  idea  as  to  the  reiative  accuracy  of  our  fully  numerical  DFT  energies  we 
note  that  Becke  and  Dickson  have  taken  L  up  to  11,  but  their  reported  energy  for,  e.g., 
CO2  (see  Table  II  in  A.D.Becke  and  R.M.Dickson,  J.Chem.Phys.,  89  (1988)  2993) 
deviates  by  0.016  hartree  from  the  exact  vaiue,  compared  to  our  error  of  0.000062 
hartree  (Tabie  6).  The  current  default  in  DMol  is  L  =  4. 


3.  Other  Projects 

In  addition  to  our  two  major  research  topics,  a  number  of  other  research  projects  in  the 

Pulay  group  have  been  supported  from  our  Air  Force  grant. 

•  We  have  implemented  a  new  Gaussian-weighted  operator,  originally  developed  by 
Rassolov  and  Chipman  for  MCSCF  wavefunctions  (V.A.Rassolov  and  D.A.Chipman, 
J.Chem.Phys.  105  (1996)  1470),  for  the  calculation  of  spin  densities  at  nuclei  in 
density  functional  theory.  Applications  to  first  row  atoms  and  some  diatomic  and 
small  polyatomic  molecules  show  good  agreement  with  experiment. 

•  We  have  carried  out  a  thorough  investigation  into  the  inner-hydrogen  migration  and 
ground-state  structure  of  porphycene  (a  structural  isomer  of  porphyrin).  We  were 
able  to  convincingly  show  via  good  quality  density  functional  calculations  and  an 
SQM  force  field  analysis  of  the  vibrational  spectrum,  that  the  ground-state  structure 
is  the  C2h  trans  isomer  and  a  previously  proposed  scheme  for  the  inner-hydrogen 
migration  involving  both  the  trans  and  cis  isomers  was  not  viable. 

•  We  have  investigated  the  structure  and  IR  spectra  of  a  number  of  metal  tris- 
acetylacetonates,  correcting  some  of  the  experimental  assignments  and  predicting 
the  geometry  and  vibrational  spectrum  of  several  species  (e.g.,  the  scandium 
compound),  which,  to  our  knowledge,  have  not  been  experimentally  determined. 
Excellent  agreement  was  obtained  with  experimental  IR  spectra  for  the  Al,  Fe  and 
Cr  compounds. 
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